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ABSTRACT 

The observables in a strong gravitational lens are usually just the image positions and 
sometimes the flux ratios. We develop a new and simple algorithm which allows a set 
of models to be fitted exactly to the observations. Taking our cue from the strong 
body of evidence that early-type galaxies are close to isothermal, we assume that the 
lens is scale-free with a flat rotation curve. External shear can be easily included. 
Our algorithm allows full flexibility as regards the angular structure of the lensing 
potential. Importantly, all the free parameters enter linearly into the model and so the 
lens and flux ratio equations can always be solved by straightforward matrix inversion. 
The models are only restricted by the fact that the surface mass density must be 
positive. 

We use this new algorithm to examine some of the claims made for anomalous 
flux ratios. It has been argued that such anomalies betray the presence of substantial 
amounts of substructure in the lensing galaxy. We demonstrate by explicit construction 
that some of the lens systems for which substructure has been claimed can be well-fit 
by smooth lens models. This is especially the case when the systematic errors in the 
flux ratios (caused by microlensing or differential extinction) are taken into account. 
However, there is certainly one system (B 1422+231) for which the existing smooth 
models are definitely inadequate and for which substructure may be implicated. 

Within a few tens of kpc of the lensing galaxy centre, dynamical friction and tidal 
disruption are known to be very efficient at dissolving any substructure. Very little 
substructure is projected within the Einstein radius. The numbers of strong lenses for 
which substructure is currently being claimed may be so large that this contradicts 
rather than supports cold dark matter theories. 

Key words: gravitational lensing - galaxies: structure - galaxies: elliptical - dark 
matter 



1 INTRODUCTION 

The fitting of models to observational data has always been a 
major concern in strong gravitational lensing (e.g., Schechter 
2000). In most cases, the lens is only composed of 2 (a "dou- 
blet") or 4 (a "quadruplet") point-like images. For the dou- 
blets, the observable constraints are the four relative coor- 
dinates of the images with respect to the lensing galaxy, 
together with the one flux ratio of the first image to the 
second. For the quadruplets, this becomes eight relative co- 
ordinates and three flux ratios. Only for a handful of lens 
systems are time delays available. It therefore happens very 
frequently that lens models have more degrees of freedom 
than the number of constraining observations. In fact, the 
situation is often even worse than we have just described. 
Sometimes the centre of the lensing galaxy itself cannot be 
reliably identified. Sometimes a lens system is complicated 



by the possible effects of nearby bright galaxies. Very often, 
the flux ratios are untrustworthy, either because of differen- 
tial extinction in the optical bands (e.g., Falco et al. 1999) 
or because of microlensing in the optical and radio (e.g., Ir- 
win et al. 1989; Koopmans & de Bruyn 2000) or because of 
scintillation, scatter-broadening and free-free absorption in 
the radio (e.g., Patnaik et al. 1992, Jones et al. 1996, Winn, 
Rusin & Kochanek 2003). It is therefore already clear that 
considerable caution must be exercised in interpreting the 
results of fits to gravitational lens systems. 

The approach of modellers has been largely two- 
pronged. First, general non-parametric methods have been 
developed (e.g., Saha & Wiliams 1997). These have proved 
powerful in exploring the range of degeneracies in the lens- 
ing mass that can give rise to a particular image configu- 
ration. For example, in the case of PG 1115+080, Saha & 
Williams present three different models obeying the lens- 
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ing constraints, but for which the lensing mass distribution 
resembles a face-on spiral, an edge-on disk or a flattened el- 
liptical galaxy respectively. Second, many simple paramet- 
ric models have been thoroughly explored, such as those 
based on elliptically stratified potentials (e.g., Witt f996; 
Witt & Mao f997) or densities (e.g., Kassiola & Kovner 
f993; Munoz, Kochanek & Keeton 200f). The advantage of 
this is that the models may already incorporate some of 
the known properties of nearby galaxies. However, this ap- 
proach may also be dangerous because simple ansatze like 
elliptical potentials or surface mass densities can introduce 
unexpected properties. These properties may be so severe 
that, for example, the flux ratios may not be well fitted and 
wrong conclusions may be drawn. This can be seen most 
clearly in the case of elliptical potentials, in which there are 
strong constraints on the flux ratios (e.g., Witt & Mao 2000; 
Hunter & Evans 2001). 

All this attains added significance in the light of recent 
claims of evidence for substructure from "anomalous flux ra- 
tios" in strong lensing. Mao & Schneider (1998) were the first 
to point out that substructure may be needed to explain the 
flux ratios in some cases. The instance that they selected, the 
quadruplet B 1422+231, has three highly magnified bright 
images and one much fainter image. Three highly magnified 
images occur generically near a cusp, and a Taylor expan- 
sion gives a universal relationship that the sum of the fluxes 
of the two outer images should equal the flux of the mid- 
dle image. This is strongly violated in B 1422+231 leaving 
substructure on top of the smooth model as the believable 
culprit. This result is supported by the detailed models of 
B 1422+231 by Bradac et al. (2002) which find that the dis- 
crepancy requires substructure on the mass scale ~ 1O 6 M . 
Recently, Dalai & Kochanek (2002) looked at seven predom- 
inantly radio quadruplets and argued that the flux ratios 
were anomalous by comparison with those expected for sim- 
ple isothermal lenses with external shear. They claimed that 
the anomalous flux ratios implied the existence of ~ 2 per 
cent of the mass of the lensing galaxy in substructure. Met- 
calf & Zhao (2002) looked at five optical quadruplets and 
similarly claimed evidence for anomalous flux ratios on com- 
parison with those expected for simple elliptical power-law 
potentials with shear. The problem with this procedure is, 
of course, that anomalous flux ratios may not be the result 
of substructure at all, but may simply reflect deficiencies in 
the modelling. 

This motivates us to introduce in sections 2 and 3 a 
new approach for fitting which can incorporate the flux ra- 
tios at outset and which can permit an arbitrary azimuthal 
dependence for the surface density. Importantly, all the free 
parameters enter linearly into the models and so the lens and 
flux ratio equations can always be solved by straightforward 
matrix inversion. The models are only restricted by the fact 
that the surface mass density must be positive. This algo- 
rithm is used in section 4 to assess the evidence for anoma- 
lous flux ratios. Readers who are mainly interested in this 
application may skip section 3 entirely on the first reading. 
This section is rather mathematical and explicitly derives 
the linear equations for the parameters to be fitted. 



2 SCALE-FREE MODELS 

2.1 Potential and Surface Mass Density 

Scale-free galaxy models with arbitrary power-law fall-off 
are widely used in galactic astronomy and dynamics (e.g., 
Toomre 1982; Evans 1994; Evans, Carollo & de Zeeuw 2000). 
The isophotes of a scale-free galaxy have the same shape at 
every radius and are completely described by a shape func- 
tion G{9), which depends only on the position angle 9 with 
respect to the major axis. In such galaxies, the convergence 
k and deflection potential <f> are proportional to a power of 
r, namely 

K =±G(0)r f '- 2 , ct> = r 9 F(e). (1) 

Here, (r, 9) are familiar polar coordinates in the lens plane. 
The power is restricted to lie in the range < (3 < 2 for 
galaxy models. For /3 < 0, the mass at the centre of the 
galaxy does not converge. For (3 > 2, the surface mass den- 
sity increases with radius, which is unrealistic. The impor- 
tant (3=1 case corresponds to an everywhere flat rotation 
curve. 

In gravitational lensing, some simple scale-free poten- 
tials were already investigated by Kassiola & Kovner (1993). 
They were followed by Witt, Mao & Keeton (2000), who 
pointed out that the time delay is independent of the angu- 
lar function F(9) for the special case (3 = 1. Many applica- 
tions to this general model followed (see Evans & Witt 2001; 
Zhao & Pronk 2001; Wucknitz 2002; Cardone et al. 2002). 

From Poisson's Law, we have 

G{9) = P 2 F(9) + F"(9). (2) 

Given F(9), it is straightforward to generate G(9). However, 
it is also possible to solve this ordinary differential equation 
for F(9) in terms of G(9), as shown in Appendix A. This is 
useful, as it gives the scale-free potential corresponding to a 
given set of scale-free equidensity contours. 

The critical curves can be given analytically. Using the 
determinant of the Jacobian, we obtain 

det J(r,9) = (f3~l)r 20 - 4 {[3FF" -((3-l)F' 2 +(3 2 F 2 ] 
- r P- 2 {F" + f3 2 F] + 1 = 0. (3) 

Setting p = r 13 ' 2 yields a quadratic equation in p which 
can easily be solved. The two roots correspond to the radial 
and the tangential critical curves respectively. They can be 
mapped onto the radial and tangential caustics with the lens 
equation using the components of the deflection angle (A2). 
For the flat rotation curve case {(3 = 1), the equation for the 
critical curve becomes linear and thus the caustic network 
is readily established, as shown in Evans & Witt (2001). 

2.2 Fourier Expansions 

Since F(9) and G(9) are periodic with period 2-k, we can 
expand them as Fourier series (see e.g. Bronshtein & Se- 
mendyayev 1998). For convenience, we start with the poten- 
tial F(9) and write 

oc 

F(9) = ^ + Y^Wk cos(fctf) + b h sin(fc0)]. (4) 

fe=i 

Using eq.(2), we obtain now for G(9) 
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Figure 1. The two plots show truncated Fourier series approximations to the equidcnsity contours of an exactly elliptical E3 (left panel) 
and an E7 (right) galaxy. For moderate flattcnings such as E3, excellent results are obtained with only the first three non- vanishing 
terms in the expansion. Even for the most highly flattened configurations such as E7, the approximation becomes very accurate when 
the first seven non-vanishing terms are included. 



G(6) 



2 

aoP 



+> M/T-fc 2 ) cos(k8)+b k (p'-k') sin(fc0)].(5) 



The coefficients in the Fourier series must obey the condition 
G{9) > for all 9 to obtain a positive surface mass density. 
After constructing a fit, this is easy to check a posteriori. 

As an example, let us consider an elliptical surface mass 
density with a flat rotation curve (/3 = 1), which has 



G{6) = A(cos 2 6» + ^ 2 sin 2 ( 



-1/2 



(6) 

Since G(9) — G(—6) is an even function, all b k must vanish 
(fefc = 0). Since G{9) = G(n — 9), all a k with k equal to an 
odd integer must also vanish. So, the only terms that remain 
in the Fourier series are the even integer cosine coefficients 
(cf. also Appendix B) . Figure 1 shows how few terms are re- 
ally required for accurate approximation. For an E3 galaxy, 
a Fourier series truncated after the first three non-vanishing 
terms already gives an excellent approximation to the true 
equidensity contours. Even for a highly flattened E7 galaxy, 
only the first seven non-vanishing terms are needed. This 
demonstrates that the Fourier expansions can give good re- 
sults with few terms. 

The advantage of the Fourier expansion is that it makes 
it very easy to solve the lens equation. For this reason, we 
list the Fourier coefficients for the popular elliptic density 
and potential models with flat rotation curves in Appendix 
B (although we do not need to use them in the main body 
of this paper). 

2.3 The Image Positions and Fluxes 

The lens equation relates the position of the source (£, rj) to 
the positions of the image (a:, y) through the derivative of 
the lensing potential (e.g., Schneider, Ehlers & Falco 1992) 



dx ' 



v = y- 



dy' 



(7) 



Using polar coordinates in the image plane, the lens equation 
can be recast as 



r cos U — r 



13-1 



rj — r sin u — r 



[/?cos0F(0)-sin0F'(0)], (8) 
\J3 sin OF (6) + cos 9F 1 (9)] , (9) 



where we can express F(6) and F'(0) by Fourier series. Let 
us assume that we have measured the image positions (re, 9i) 
(in polar coordinates) of a gravitational lens and inserted 
them into the lens equation. Immediately, we notice that all 
unknown quantities (i.e. £, rj, and all a,i and bi) enter linearly 
into the lens equation, except for j3. 

Let us now add in the constraint that the magnification 
ratios of the images are equal to the observed values. Sup- 
pose that the ratio of the magnification of the £th image to 
the fcth is measured to be f k e, then 



fkt 



det J(r e ,8 e 
det J(r k ,9 k 



(10) 



where det J is given in eq. (3). We observe that the coeffi- 
cients cii and bi enter in simple mixed quadratic form into 
the equations. However, for the special case of a flat rotation 
curve (fi) = 1), all unknown quantities enter linearly into the 
equation for the magnification relation. 

In other words, in the astrophysically important flat 
rotation curve case, the Fourier coefficients describing the 
lensing galaxy are related to the observables by a simple ma- 
trix equation. This makes the problem of fitting to both the 
image positions and the flux ratios a straightforward matter 
of matrix inversion. 

A number of recent investigations have exploited the 
well-known linearity of the lens equation to expand the po- 
tential in basis functions (e.g., Keeton 2001). Such methods 
can be a powerful way of exploring a large set of models, al- 
though the drawback is that some of the parameter space is 
often associated with the creation of additional images. The 
need to check for this can undermine the utility of starting 
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with such an algorithm. It is worthwhile emphasising how 
our work differs from such expansions. First, the fundamen- 
tal and new contribution of our paper is that the flux ratio 
equation is linear for a galaxy with a flat rotation curve. 
Therefore, in our algorithm both the lens equations and the 
flux ratio equations are simultaneously solved. Second, the 
caustic network is straightforwardly available using the work 
in Section 2.1, and so it is very easy to check for the creation 
of unwanted additional images. 



3 FITTING GRAVITATIONAL LENSES 

3.1 A Singular Value Decomposition Problem 

First, we insert the Fourier series of F(9) into eqs. (8)-(9) 
to establish the linear equations for fitting for /? = 1. Intro- 
ducing a useful notation, we write for the lens equation 

oo 

£M) = rms 9-^-a (9)-^[a k a k (9)+b k p k (9)l (11) 



r,(r,9) = rsme-^-a o (9)-^2[a k a k (0)+b k p k (e)l (12) 

fc=i 

where a k and b k are the Fourier coefficients of the unknown 
function F(9), which describes the angular part of the po- 
tential. For the coefficients a k and p k , we find 



a k (9) — cos9cos(k9) + k sin 9 sin(k9) , 

f3 k {0) = cos(9sin(fc(9) - k sin 8 cos(k9) , 

and for & k and p k , we find 

& k (9) — sin 6 cos(k9) — k cos 6 sin(k9), 

p k (9) = sin9sin(k9) + k cos 9 cos(k9) , 



(13) 



(14) 



for all k > 0. The lens equations (8)-(9) must hold at the po- 
sitions of the n images, given in polar coordinates as (re, 9e) 
with I = l,...,n. This provides 2n constraints on the un- 
known source position (£, rj) and the unknown Fourier coef- 
ficients. 

Now let us consider the magnification of the £th image. 
For the flat rotation curve case (P = 1), we obtain 



det J = 



1 



G(0 t ) 



(15) 



pe^e re 

where pe is the parity of the image I. It is only the flux 
ratios that can be related to the observational data, not 
the magnifications themselves. For the flux ratio of image i 
compared to image k, we obtain 

- G(9e)/re 



f, 



p k ^ k 



(16) 



pefie l-G(9 k )/r k ' 

The flux ratio may be positive or negative depending on the 
combination of parities of the images. For flux ratios, we 
always compare to a reference image, which without loss of 
generality we take to be the first image. As we take ratios, 
we only need to identify images of the same parity, which in 
a quadruplet lens are usually on opposite sides of the lensing 
galaxy. Clearing the fractions, we can write the equation in 
a similar form to the lens equation 

oo 

(fu-l)rm = ^lo{9e)+Yyklk(9e)+b k 5 k {9 t )], (17) 



with 

l k {9e) 
S k (9 e ) 



(l-k 2 )[f 1 ere cos(k9 1 )-r 1 cos(fc^)] 
(l-k 2 )[f 1 ere sin(k9 1 )-r 1 sin(fc^)]. 



(18) 



This provides a further n — 1 linear constraints on the un- 
knowns. 

Now let us introduce a compact notation to facilitate 
matters. We set a k e = a k (9e), where the first index de- 
notes the Fourier component and the second index the im- 
age number. Exactly similar definitions can be made for 
Pk£,& k e, P k e,lke and 5 k e- We are now ready to set up the 
matrix equation for the fitting to the positions of the im- 
ages and the flux ratios. Let us write 



Cx = d 

where 

d = (re cos ( 



.,r e siaOt, (fu - l)rin, •••) 



(19) 



(20) 



is a vector of 3n— 1 observables. Here, I runs from 1 to n 
for a lens system with n images. There is always one less 
equation for the flux constraints than for the image position 
constraints. The vector 



x = (£,77,00/2, 01,61,02,62, ■■■■) 



(21) 



contains the unknown quantities, namely the source coordi- 
nates (£,f)), and the Fourier coefficients. There are in prin- 
ciple infinitely many such coefficients, but in practice the 
Fourier series is usually terminated so that a; is a vector 
of 3n — 1 components as well. Finally, for the matrix C, we 
obtain 

/ 1 aoe ecu Pie ciu P21 



C = 



1 &oe &u P\ 



7( W 7i£ 8_ 



e &2e Pie 
e J2e &2e 



(22) 



Provided we can measure all n image positions and n— 1 flux 
ratios, then C has 3n — 1 rows. The number of columns is 
in principle arbitrarily large, although in practice it usually 
makes sense to terminate the Fourier series so that C is a 
3n — 1 x 3n— 1 matrix. 

As it stands, the matrix C is always singular (det C = 
0). There exists a null space for the matrix C, i.e., there 
exists at least one non-vanishing vector xo which satisfies 
the equation Cxo = 0. In fact, the matrix C has (at least) 
a two-dimensional null space. Fortunately, the null subspacc 
can be easily constructed. Since au = 1, pu = 0, &u = 0, 
Pie = 1, 1u = and 5u = for all t, it is easy to verify 
that the null vector must be of the form xq = A1V1 + A2i»2 
with vi = (1, 0, 0, -1, 0, ...) T and v 2 = (0, 1, 0, 0, -1, 0, ...) T . 

From the physical point of view, the Fourier coefficients 
ai and 61 simply produce a shift or transformation of the 
source position (£, rf). However, both coefficients do not con- 
tribute to the surface mass density for the special case of 
p = 1 since ai cos 9 and 61 sin 9 are solutions of the homoge- 
nous differential equation F(9) + F"(9) = (cf. eq.(2)). We 
can excise this degeneracy by simply setting ai = 61 = 0, 
which corresponds to the choice of origin. Now we can just 
remove the 4th and 5th column of the matrix C and add two 
higher degree Fourier coefficients to the problem to main- 
tain C as a 3n— 1 x 3n— 1 matrix. The new matrix is usually 
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—Aa 


A<5 


Radio fluxes 


Mid-infrared flux 


H band 


I band 




(in ") 


(in ") 


(in fj, Jy) 


fractions 


(in mags) 


(in mags) 


A 


0.0 


0.0 


65.5 ±8.4 


0.27 ±0.02 


14.96 ±0.06 


15.92 ±0.12 


B 


0.673 ± 0.003 


1.697 ±0.003 


64.2 ± 8.4 


0.30 ±0.02 


15.46 ±0.02 


17.21 ±0.11 


C 


-0.635 ±0.003 


1.210 ±0.003 


26.5 ± 8.4 


0.16 ±0.02 


15.71 ±0.03 


16.77 ±0.12 


D 


0.866 ± 0.003 


0.528 ±0.003 


59.4 ±8.4 


0.27 ±0.02 


16.00 ±0.04 


17.39 ±0.04 


G 


0.075 ± 0.004 


0.939 ±0.003 











Table 1. Observational data on the Einstein Cross. The optical positions are taken from the CASTLES survey, as are the H and I band 
fluxes. The radio fluxes are provided by Falco ct al. (1996). The mid-infrared flux fractions arc given in Agol, Jones & Blacs (2000). 



Band 




V 




a 2 


f>2 


CL3 


f>3 


a 4 




a 5 


f>5 


Radio 


0.0698 


-0.0134 


1.7748 


-0.0422 


0.0430 


0.0004 


-0.0017 


0.0007 


0.0012 


0.0000 


0.0008 


Mid-IR 


0.0662 


-0.0131 


1.7711 


-0.0417 


0.0417 


0.0000 


0.0002 


0.0009 


0.0005 


-0.0002 


0.0001 


I 


0.0689 


-0.0142 


1.7670 


-0.0368 


0.0477 


-0.0043 


0.0018 


-0.0016 


-0.0018 


0.0029 


0.0029 


H 


0.0640 


-0.0132 


1.7718 


-0.0392 


0.0448 


-0.0013 


0.0011 


0.0002 


0.0004 


0.0019 


0.0025 



Table 2. Fourier coefficients of the solutions for the Einstein Cross. They correspond to the bold lines in Figure 2, for which the positions 
and flux ratios are reproduced within the errors and are given in the same coordinate system as the data in Table 1 but translated to 
the lens centre. 




Figure 2. Each panel shows in a dashed line the critical curve (which is also an equidensity contour when the rotation curve is flat) of a 
lens model for the Einstein Cross that exactly fits the image positions and the flux ratios. The full curve shows the equidensity contour 
of a model that fits the data within the uncertainties. Working clockwise from the top left, the panels use the flux ratios in the radio, 
mid infrared, H band and I band respectively. The locations of the four images and the lensing galaxy are also marked. 
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-Act 


fin "1 


AS (ir 


l 


Radio Fluxes 




Predicted 


Observed 


Predicted 


Observed 


Predicted 


Observed 


A 


0.0 


0.0 


0.0 


0.0 


65.5 


65.5 


B 


0.673 


0.673 


1.697 


1.697 


61.1 


64.2 


C 


-0.635 


-0.635 


1.210 


1.210 


27.9 


26.5 


D 


0.866 


0.866 


0.528 


0.528 


56.6 


59.4 


G 


0.070 


(0.075) 


0.939 


(0.939) 








-Aa (in ") 


AS (ir 


") 


Mid-infrared flux fractions 




Predicted 


Observed 


Predicted 


Observed 


Predicted 


Observed 


A 


0.0 


0.0 


0.0 


0.0 


0.27 


0.27 


B 


0.673 


0.673 


1.697 


1.697 


0.27 


0.30 


C 


-0.635 


-0.635 


1.210 


1.210 


0.15 


0.16 


D 


0.866 


0.866 


0.528 


0.528 


0.30 


0.27 


G 


0.075 


(0.075) 


0.939 


(0.939) 







Table 3. Predicted and observed quantities for two of the solutions of the Einstein Cross given in Figure 2. The upper (or lower ) 
panel corresponds to the lens model with equidensity contour given by the bold curve in the panel for the radio (or mid-infrared fluxes). 
The positions —Aa and A<5 are computed relative to image A. Notice that all the data are well-reproduced within even the formal 
uncertainties and the flux ratios arc not anomalous. 





Figure 3. The panels show the caustics of the Einstein Cross corresponding to the exact solutions (dashed curves in Fig. 2). The 
computed position of the lensed quasar is indicated. It always lies in the region of 4 images. For the I and H band, there are butterfly 
and swallowtail cusps (small regions of the source plane where 6 images may occur). 
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non-singular and so a Gauss- Jordan elimination can formally 
solve the problem. However, it is better to use singular value 
decomposition (SVD) to solve for the unknown vector x (cf. 
Press et al. 1999). In this way, we are always guaranteed to 
find a numerically stable solution. 

The singular value decomposition of the matrix C is 

C = U ■ W ■ V T , (23) 

where the 3n— 1 x 3n— 1 matrices U and V are orthonormal, 
and the matrix W is a diagonal matrix of singular values 
(wi, W2, ■ W3n-i)- As usual in a singular value decompo- 
sition (SVD), any very small singular values are removed 
(see Press et al. 1989) and this must be accompanied by a 
corresponding reduction in the number of unknown Fourier 
coefficients. Typically, any singular value whose ratio to the 
largest singular value is less than 1CF 4 is removed. After 
deletion, let us suppose that m singular values remain, where 
m < 3n — 1. If m = 3n — 1, then the SVD solution is exact; 
if singular values are removed, then the solution is usually 
very good. 

Let us remark here that one can add further higher 
degree Fourier coefficients and then the fitting problem be- 
comes underdetermined. We can obtain a whole space of 
models all of which would equally well reproduce the ob- 
served data. It is clear - even with the severe restriction 
that we have made to a scale-free model with a flat rota- 
tion curve - that the variety of models that fit the data is 
extremely large. This cautions us against making rash state- 
ments based on goodness-of-fit to a single model. 

3.2 Uncertainties 

It is of modest interest to provide a model that fits the 
observational data exactly, as the uncertainties in the flux 
ratios may be substantial. It is of much greater interest to 
be able to find a set of models that can reproduce the data 
within the errors. 

The relative image positions are usually known to ex- 
tremely high accuracy. For lenses observed with the Hubble 
Space Telescope, the error in the relative astrometry is only 
~ 3 mas. From the point of view of modelling, however, the 
important quantities are the relative positions of the images 
with respect to the centre of the mass distribution of the 
lensing galaxy. These coordinates are less accurately deter- 
mined. The centre of the light distribution of the lensing 
galaxy is sometimes poorly known (especially if it is at high 
redshift). In any case, there is no guarantee that the center 
of the light corresponds to the centre of the mass. In our 
modelling, we always reproduce the relative image positions 
exactly, but we allow for an error of no more than 50 mas 
in the position of the center of the lens. 

For the flux ratios, we distinguish between the radio 
and the optical data. The flux measurements in the radio 
may be affected by scintillation or scatter-broadening or 
free-free absorption. An isolated example of radio microlens- 
ing is known (Koopmans & de Bruyn 2000), but the effects 
of microlensing are not generally believed to be substantial 
in the radio. Whenever we use radio data, we assume that 
the typical uncertainty in the flux ratio is no more than 5 
per cent. However, there are much greater problems in the 
optical bands. Each image might be affected by differen- 
tial extinction or by microlensing. Microlensing can cause 



an image to be below or above its average magnification for 
decades (see Witt & Mao 1994 for an extensive discussion 
of the effects of microlensing). The recent years have seen 
monitoring campaigns of both Q 2237+030 and Q 0957+561, 
which provide irrefutable evidence for microlensing (Col- 
ley et al. 2000; Wozniak et al. 2002; Schmidt et al. 2002). 
There are at least a further five multiply imaged quasars for 
which evidence exists implicating microlensing in the opti- 
cal (HE 1104-1805, PG 1115+080, H 1413+117, B 0218+357, 
B 1600+434). Therefore, the optical flux ratios are more un- 
certain and we assume that they have an error of < 10 per 
cent. 



3.3 Position Angle of the Lensing Galaxy Axes 

The axes of the local right ascension and declination coor- 
dinates may be rotated by an arbitrary angle with respect 
to the major and minor axes of the lensing galaxy. This is 
implicitly contained in the derived Fourier coefficients. 
We can expand the Fourier series as 

oo 

^^[cifc cos(fc6<) + bk sin(fc#)] 
fe=i 

oo 

= ^2[akCO8(k6 + k6 R )+bksm(k0 + k0R)], (24) 
fc=i 

where 9 R is the rotation angle and a k and b k are the Fourier 
coefficients in the rotated coordinate system. It is easy to 
verify that 

o-k = a k cos k9 R + b k sin k9 R 

b k = — a fe sin k9 R + b k cos k9 R . (25) 

If the lensing galaxy has a major axis, the following condi- 
tion must be satisfied: 

G'(0 R ) «O«G'(0fl + 18O o ). (26) 

If the lensing galaxy looks regular, it probably has, in ad- 
dition to a major axis, an approximate reflection symmetry 
G(9) « G(-9), or equivalents F{6) « F(-0) . Therefore, 
in the rotated system, we expect b k « for all k. In other 
words, the rotated system is defined by 

n 

^~^(bfc ) 2 = minimum. (27) 
fe=i 

Solving this equation yields the rotation angle of the major 
axis of the lensing galaxy 

n n 

2^ka k b k cos(2k9 R ) =^k{a k - b 2 k )sm(2k9 R ). (28) 

k=l k=l 

This equation must in general be solved numerically. For 
galaxy-like mass distributions, the dominant higher-order 
terms in the Fourier expansion are those with k = 2 and 
then the solution is analytic, namely 

tan20 H =— . (29) 

For applications to individual lenses, we note that the posi- 
tion angle as conventionally defined by astronomers (North 
through East) is 90° - 6 R . 
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—Act 


AS 


I band (K93) 


H band 


V band 


I band 




(in ") 


(in ") 


(in mags) 


(in mags) 


(in mags) 


(in mags) 


Al 


-1.328 ±0.004 


-2.034 ± 0.004 


16.12 


15.71 ± 0.03 


16.90 ±0.11 


16.42 ± 0.02 


A2 


-1.477 ±0.004 


-1.576 ±0.003 


16.51 


16.21 ± 0.04 


17.62 ± 0.09 


16.85 ± 0.03 


B 


0.341 ±0.003 


-1.961 ±0.004 


18.08 


17.70 ± 0.05 


18.39 ± 0.45 


17.91 ± 0.39 


C 








17.58 


17.23 ± 0.04 


18.95 ± 0.32 


18.37 ±0.34 


G 


-0.381 ±0.003 


-1.344 ±0.003 











Table 4. Observational data on the quadruplet PG1115±080. The optical positions are taken from the CASTLES survey, as arc the 
V,H and I band fluxes. For comparison, we also give the earlier I band fluxes found by Kristian ct al. (1993). 



Band 


s 


V 


a 


a 2 


£>2 


6E3 


h 


a 4 


64 


H 


-0.0273 


0.1122 


2.2946 


-0.0029 


0.0052 


-0.0021 


0.0027 


-0.0007 


-0.0002 


I (K93) 


-0.0277 


0.1114 


2.3039 


0.0007 


0.0059 


0.0013 


0.0060 


0.0000 


0.0020 


V 


-0.0238 


0.1175 


2.2871 


-0.0081 


0.0035 


-0.0054 


-0.0001 


-0.0022 


-0.0020 


I 


-0.0218 


0.1145 


2.2880 


-0.0064 


0.0034 


-0.0049 


0.0001 


-0.0016 


-0.0012 



Table 5. Fourier coefficients of the solutions for PG 1115+080, assuming a constant external shear with a magnitude of 0.103 and a 
direction of 65.8° measured from west to north (following Schcchtcr ct al. 1997). The coefficients are computed in the same coordinate 
system as the data in Table 4 but translated to the lens centre. 



3.4 Mass inside the Einstein Ring 

Evans & Witt (2001) showed that the mass inside the Ein- 
stein ring must be very close to the mass inside the critical 
curve for the special case of (3 = 1, and can be written as 

Me w M crit . curvc = -L / G 2 (8)d6. (30) 
27V Jo 

We therefore deduce that 

M ^T + X> 2 - 1)2(ai + b|) - (31) 

To apply this to a real lens, we have to scale the result by 
the factor E cr it-Dd containing the normalized surface mass 
density which is is given by 

c 2 D B 

Ecrit = i-kGdId^ (32) 

Here, D^,D S and D^ s are the distance to the deflector, the 
distance to the source and the distance from the deflector 
to the source respectively. Given the redshift of the lens and 
the source, this (32) enables the projected mass within the 
Einstein ring to be estimated. 

4 AN APPLICATION: ANOMALOUS FLUX 
RATIOS 

In this section, we consider three of the lens systems for 
which substructure has been claimed either by Metcalf & 
Zhao (2002) or by Dalai & Kochanek (2002) or by Chiba 
(2002). 

4.1 Q 2237+030 (The Einstein Cross) 

The data on the four images of Q 2237±030 (the Einstein 
Cross) are listed in Table 1. Much of this comes from the 



CASTLES survey *. Q2237±030 is an unusual lens because 
so much more information is available for the lensing galaxy 
than is customary. In particular, the redshift of the lens is 
so small (z = 0.039) that the galaxy's light distribution can 
be measured (e.g., Wyithe et al. 2002). The lens is a face-on 
barred galaxy (e.g., Schmidt, Webster & Lewis 1998). The 
optical fluxes are known to be affected by microlensing (Ir- 
win et al. 1989; Wambsganss, Schneider & Paczyhski 1990; 
Wozniak et al. 2000) , which causes each image to differ from 
the average magnification. Typically, we expect that the im- 
ages are very likely below the average magnification if they 
are in a quiescent state (for years), while the images are very 
likely above the average magnification if they are in a more 
active state (cf. Witt & Mao 1994). 

It is often argued that the radio and mid-infrared fluxes 
are more reliable because the source emitting region is more 
extended and therefore less affected by microlensing. If this 
is the case, then gravitational lens models should use the 
radio and mid-infrared flux ratios in preference to optical. 
The flux ratios in the radio were measured by Falco et al. 
(1996) with a signal-to- noise of ~ 2 — 4. The mid-infrared 
flux ratios have recently become available thanks to Agol, 
Jones & Blaes (2000). The flux ratios in radio and mid- 
infrared are in good agreement, but differ from the optical 
flux ratios. 

Figure 2 shows the critical curves for four models of 
Q 2237+030. From eq. (15), the critical curve is also an 
equidensity contour when the rotation curve is flat. As the 
models are scale-free, the equidensity contours retain the 
same shape independent of position, although the size of 
course varies. In each case, the dashed curve describes a 
model that exactly reproduces the image positions and flux 
ratios. It is found by solution of the matrix equation (19) 
using singular value decomposition (SVD). Figure 3 shows 

* http://cfa-www.harvard.edu/castlcs/ 
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— Aa 


fin "1 


A<5 (in ") 


H band fluxes (in mags) 




Predicted 


Observed 


Predicted 


Observed 


Predicted 


Observed 


Al 


-1.324 


-1.328 


-2.035 


-2.034 


15.71 


15.71 


A2 


-1.481 


-1.477 


-1.575 


-1.576 


16.10 


16.21 


B 


0.341 


0.341 


-1.961 


-1.961 


17.59 


17.70 


C 














17.33 


17.23 


G 


-0.391 


(-0.381) 


-1.334 


(-1.344) 








-Aa (in ") 


AS (in ") 


I band fluxes of K93 (in mags) 




Predicted 


Observed 


Predicted 


Observed 


Predicted 


Observed 


Al 


-1.324 


-1.328 


-2.035 


-2.034 


16.12 


16.12 


A2 


-1.481 


-1.477 


-1.575 


-1.576 


16.40 


16.51 


B 


0.341 


0.341 


-1.961 


-1.961 


18.18 


18.08 


C 














17.47 


17.58 


G 


-0.386 


(-0.381) 


-1.329 


(-1.344) 







Table 6. Predicted and observed quantities for two of the solutions given in Figure 4. The upper (or lower) panel corresponds to the 
lens model with equidensity contour given by the bold curve in the panel for the H band (or I band of K93) fluxes. The positions — Aa 
and AS arc computed relative to image C. Notice that astrometric data are well-reproduced within even the formal uncertainties whilst 
the flux ratios are not anomalous, once an uncertainty of 10 per cent is allowed for the effects of microlensing. 




Figure 4. The panels show the equidensity contours (re = 1/2) for lens models of PG 1115+080. The bold curves reproduce the relative 
astrometry within the errors and the measured optical flux ratios to better than 10 per cent. As a constant external shear is assumed, 
the equidensity contour is no longer a critical curve. The locations of the four images and the lensing galaxy are also marked. 
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Figure 5. The panels show the caustics of the lens model for PG 1115+080 corresponding to the solutions in Fig. 4. The computed 
position of the lensed quasar is indicated. It always lies in the region of 4 images. There are no regions in which higher order imaging 
occurs (butterfly and swallowtail cusps). 



the corresponding caustics of the exact solutions of the crit- 
ical curves (dashed lines) . The more irregular critical curves 
(I and H band data) correspond to caustics with over-folds 
like swallowtails and butterflies. However, in each case, the 
source position still lies inside the four image region. Un- 
wanted extra images are not produced. 

We assume that the measured flux ratios in the optical 
band are in error by < 10 per cent due to microlensing. The 
relative positions of the image are always maintained within 
the formal errors (given in Table 1), but the center of the 
lensing galaxy may be displaced by up to 50 mas. We find 
the model for which 

bi + J2^ + b ^ ( 33 ) 

i>3 

is minimised. This suppresses Fourier components higher 
than the bar-like m = 2 component. In other words, out 
of all the solutions that reproduce the data, we choose the 
one that most looks like a galaxy. The models so produced 
are shown in bold lines in Figure 2; the details of the solu- 
tions are listed in Table 2. Finally, for two of the bold line 
solutions, we show how the image positions and the flux ra- 
tios are reproduced in Table 3. All the observable quantities 
are extremely well reproduced by the model, including the 
flux ratios. 



Two points emerge clearly from this. First, the equiden- 
sity contours implied by reproducing the exact radio or mid- 
infrared fluxes are already in good agreement. The position 
angle of the lensing galaxy as inferred from eq. (26) using 
the radio or mid-infrared solutions is ~ 67.3° (measured 
from North through East). For comparison, Trott & Web- 
ster (2002) report that the position angle of the bar in the 
lensing galaxy is ~ 39°, while the major axis of the galaxy 
in its outer parts is ~ 77°. Second, the equidensity con- 
tours implied by exactly reproducing the I and the H band 
optical flux ratios are distorted. However, there are more 
regular solutions once an allowance is made for the effects 
of microlensing on the flux ratios. 

Distortions of isophotes from elliptical form are gen- 
erally given in terms of Bender et al.'s (1989) coefficients. 
These are usually calculated for early-type galaxies and have 
absolute values less than 0.05 (see e.g., Binney & Merrificld 
1998, section 4.3). In fact, barred galaxies, such as the lens 
in Q 2237+030, may have substantially larger deviations, 
which moreover depend on position along the major axis. 
The method of calculating these coefficients for densities 
given as Fourier expansions is elaborated in Appendix C; we 
quote the results here. For the radio contours, Bender's co- 
efficients are af /a = —0.008 and af /oo = —0.005; for the 
mid-infrared, they are af /oo = 0.001 and af /do = —0.011. 
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These values are well within the observed range, empha- 
sising that our solutions for the shapes of the equiden- 
sity contours are realistic. By comparison, for the distorted 
I band contours, the coefficients are af/ao = 0.012 and 
of /a = -0.034. 

All this provides strong evidence for the point of view 
advanced in Agol et al. (2000) that the radio and mid- 
infrared fluxes are trustworthy and that the optical flux 
ratios are discrepant because of microlensing. All the data 
are internally consistent and once reasonable error bars are 
placed on the flux ratios, there is no need to invoke sub- 
structure to explain the data (c.f. Metcalf & Zhao 2002). 

4.2 PG 1115+080 

Let us now examine the lens system PG 1115+080, which 
has been claimed to provide evidence for substructure by 
Dalai & Kochanek (2002), Metcalf & Zhao (2002) and Chiba 
(2002). The observational data on PG 1115+080 are listed 
in Table 4. The I, V and H band fluxes are provided from the 
CASTLES survey. There is also earlier data on the I band 
fluxes thanks to Kristian et al. (1993, hereafter K93). 

There is already evidence that the flux ratios may come 
with substantial uncertainties. For example, Foy, Bonneau 
& Blazit (1985) found component A2 of PG 1115+080 was 
brighter than Al in 1984 March-May, whereas most subse- 
quent investigators have found A2 to be fainter than Al. 
The conclusion advanced by Christian, Crabtree & Waddell 
(1987) is that the images in PG 1115+080 are being affected 
either by microlensing or by time-delayed intrinsic variations 
in the quasar or both. Witt, Mao & Schechter (1995) showed 
that the likelihood of observing microlensing in images Al 
and A2 is much higher than in C and D. Notice too that 
the CASTLES and the Kristian et al (1993) I band fluxes 
in Table 4 are not consistent within the stated error bars, 
the discrepancies almost certainly being due to microlensing. 
Radio emission would be less affected by microlensing, but 
PG 1115+080 is radio-quiet and so the optical fluxes have 
to be used. 

Impey et al. (1998) have carried out deep Hubble Space 
Telescope imaging of PG 1115+080 in the infrared. They 
found that the primary lens is a nearly round elliptical 
galaxy with an ellipticity of less than 0.07. The lens is part 
of a group of galaxies, whose effects must be included in 
successful models. The group is modelled by a uniform ex- 
ternal shear with a magnitude of 0.103 and a direction of 
65.8° measured from west to north, following Schechter et 
al. (1997). The extension of our SVD method to the case 
of uniform shear is sketched in Appendix D. We implement 
our fitting procedure, deleting two singular values for numer- 
ical stability and finding the solution that minimises (33), 
but still reproduces the data within our assumed errors. The 
panels of Figure 4 show the equidensity contous of the result- 
ing models. Again, the relative astrometry of the images is 
recovered to within the formal error (~ 4 mas), the lens cen- 
tre may be misligned by up to 50 mas, while the optical flux 
ratios are reproduced to within 10 per cent. The comparison 
between model predictions and observed quantities for the 
H and 1 band (K93) solutions are given in Table 6. The con- 
tours implied by reproducing the H or I band (K93) fluxes 
are in good agreement and look plausibly like the projected 
mass distribution of a very round early-type galaxy, in ac- 



cordance with Impey et al.'s (1998) imaging. There is some 
slight distortion in the equidensity contours that reproduce 
the V or I band (CASTLES) flux ratios. The caustic net- 
works are shown in Figure 5 to demonstrate that unwanted 
extra images are not formed. 

Again, the shapes of the equidensity contours in Fig- 
ure 4 do exhibit small deviations from pure elliptical form, 
which can be quantified by Bender et al.'s (1989) coefficients. 
Using the method of Appendix C, we find af/ao = 0.009 
and af/ao = -0.004 for the H band and af/ao = -0.015 
and af/ao = —0.004 for the I band (K93) contours. There 
is some data on the shapes of dark haloes of spiral galaxies 
from numerical simulations (e.g., Heyl, Hernquist & Spergel 
1994), which supports the notion that 0.3 1 < 0.05 and 
|«4 1 < 0.05. So, the amplitudes of the higher order terms 
in our solutions for the lens galaxy do lie easily within the 
range expected for dark haloes of spiral galaxies. 

It is clear from the smooth models that we have con- 
structed that there is no really compelling evidence for 
substructure in this lens at the moment, despite abundant 
claims to the contrary (Dalai & Kochanek 2002; Metcalf & 
Zhao 2002; Chiba 2002; Kochanek & Dalai 2003). There are 
smooth and realistic models that do fit the data within the 
likely uncertainties. 



4.3 B 1422+231 

Patnaik et al. (1992) discovered the gravitational lens 
B 1422+231 in a radio survey. Table 7 lists some of the avail- 
able data on the images, largely provided by the CASTLES 
group. It was realised by both Hogg & Blandford (1994) and 
Kormann, Schneider & Bartelmann (1994) that it is difficult 
to reproduce the flux ratios of this lens system with a sim- 
ple model. Hogg & Blandford used a mixture of isothermal 
spheres and point masses whereas Kormann et al. used a 
model with external shear. The latter strategy has been fol- 
lowed by most modern modellers. 

Using our SVD method, it was impossible to reproduce 
the reported errors in the relative astrometry, together with 
the assumed uncertainties in the flux ratios (5 per cent for 
the radio, 10 per cent for the optical) at least with physi- 
cal models. Unphysical models are readily identified as the 
equidensity contours are self-intersecting. Accordingly, for 
this lens only, the errors in the relative astrometry were in- 
flated by a factor of 10 to ~ 30 mas. Figure 6 now shows 
the results of our fitting procedure in four flux bands, but 
none of the results can be said to be at all satisfactory. The 
equidensity contours are strongly distorted from those ex- 
pected in galaxy-like models. Just as alarming as the shape 
of the contours is the fact that the fitted parameters do not 
have consistent values from one band to the next. 

The conclusion to be drawn by the failure of the fitting 
procedure is that - even taking into account observational 
uncertainties - this lens system requires additional ingredi- 
ents. Possibly, as Mao & Schneider (1998) were the first to 
suggest, substructure may be causing a dramatic effect on 
the flux ratios. Possibly, as Kundic et al. (1997) suggest, the 
effects of external shear caused by nearby galaxies and clus- 
ters of galaxies may also need to be included in the analysis. 
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V band 
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(in mags) 


(in mags) 


(in mags) 


A 


0.0 


0.0 


216 


16.77 


16.43 ±0.11 


14.41 ±0.03 


B 


0.385 ± 0.003 


-0.317 ±0.003 


221 


16.45 


16.45 ±0.10 


14.29 ±0.03 


C 


0.722 ±0.003 


-1.068 ±0.003 


115 


17.25 


17.09 ±0.07 


14.98 ±0.03 


D 


-0.562 ± 0.004 


-1.120 ±0.003 


4.5 


20.40 


20.44 ± 0.06 


18.14 ±0.04 


G 


-0.375 ± 0.004 


-0.973 ± 0.004 











Table 7. Observational data on the quadruplet B 1422±231. The optical positions are taken from the CASTLES survey, as arc the H 
and V band fluxes. The radio fluxes are provided by Patnaik et al. (1992). 




Figure 6. The panels show the equidensity contours (or critical curves) for the quadruplet system B 1422±231. Working clockwise from 
the top left, the panels use the flux ratios in the H band, r band, V band and radio. This is a lens for which substructure may well be 
present, as none of our models arc at all satisfactory. 



5 DISCUSSION 

The main aim of this paper is to present a new, simple ap- 
proach to fitting a model to the data on a gravitational lens. 
Given the image positions and the flux ratios (if available), 
a model that exactly fits that data can be constructed by 
a simple matrix inversion (preferably by singular value de- 
composition, available in standard numerical libraries). All 
the fitting parameters enter linearly into the equation. No 
complicated non-linear x 2 -fitting needs to be done. 

The lensing galaxy is always assumed to be scale-free 
with a flat rotation curve. There is already a substantial 
body of evidence from fitting of lenses that early-type galax- 



ies have nearly isothermal profiles (Kochanek 1995; Munoz 
et al. 2001; Cohn et al. 2001; Rusin et al. 2002). In particular, 
the absence of central images strongly constrains the lensing 
potential of early-type galaxies to be isothermal or steeper 
and the size of any core region to be small (e.g., Rusin & Ma 
2001; Evans & Hunter 2002). There is both stellar dynam- 
ical (e.g., Gerhard et al. 2001) and X-ray (Fabbiano 1997) 
evidence that early-type galaxies have flat rotation curves 
out to at least 4 effective radii. Hence, our assumption of a 
flat rotation curve seems exactly what the data require. 

Unusually, our algorithm allows full flexibility as re- 
gards the angular structure of the lensing potential. Earlier 
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fitting procedures have allowed the lensing potential to have 
different radial structure (for example, different power-law 
profiles), but usually only simple forms of angular struc- 
ture (for example, constant external shear). In our models, 
the radial structure is always fixed as isothermal, but the 
shape of the density contours is completely arbitrary. The 
lens model has a flexible number of free parameters, and 
of course includes the isothermal elliptic density and po- 
tential models which have a distinguished history in grav- 
itational lens modelling. The fit allows direct deduction of 
the mass contour of the lensing galaxy. Although we have 
presented examples based on quadruplet lens systems, the 
algorithm can be adapted for a two image or a higher image 
model. A higher image model might be needed, for exam- 
ple, in Q 0957+561, where 10 sub-image positions can be 
detected in the radio band. Another good application might 
be MG 0414+0534, where some substructure of the images 
are detected (cf. Trotter, Winn & Hewitt 2000 and their 
fitting approach). Both lenses seem heavily distorted by ex- 
ternal shear. A shear term can be added to our algorithm. 
If the shear is known, then the lens and flux ratio equations 
remain linear. 

The flux ratios may enter directly into the fit depending 
on whether they are available or trustworthy. We stress that 
flux ratios often do not provide good model constraints. The 
reason for this is easy to understand. The flux ratios depend 
on the second derivatives of the lensing potential, whereas 
time delays and image positions are proportional to the po- 
tential and its first derivative respectively. Therefore, flux 
ratios are particularly sensitive to graininess in the gravita- 
tional potential. This often manifests itself as microlensing. 
The flux ratios are also affected by the uncertain differential 
extinction corrections that must be applied to each image. 
Radio and mid-infrared fluxes are generally more reliable 
than optical. However, effects such as scintillation, scatter- 
broadening and free-free absorption may affect the radio 
fluxes in some of the lenses. If flux ratios are incorporated 
into a fit, it is vital that there is a realistic treatment of the 
errors. 

Our fitting algorithm permits the construction of a 
whole class of degenerate models all of which can fit the im- 
age position and flux ratios exactly. In addition, the models 
are also degenerate in the time delay since the delay depends 
only on the distance of the image positions to the centre of 
the lensing galaxy (cf. Witt, Mao & Keeton 2000). Such de- 
generacy is easy to understand pictorially. Given a Fermat 
time delay surface, we need only keep the values, the deriva- 
tives and the second derivatives of the surface at the image 
positions fixed. We then have enormous freedom to move the 
surface (subject only to the constraints that no additional 
images are introduced and that no negative mass density is 
produced). 

We have used our new fitting algorithm to examine criti- 
cally some of the claims made recently for anomalous flux ra- 
tios. The procedure used by both Dalai & Kochanek (2002) 
and Metcalf & Zhao (2002) was to show that a family of 
simple models did not reproduce the flux ratios. However, 
it does not follow from this that substructure is necessary; 
it may just be that the simple model did not have enough 
flexibility to provide a good match. For example, of the five 
lenses studied by Metcalf & Zhao (2002), they themselves 
concluded that one (MG 0414+0534) could be satisfacto- 



rily explained by a smooth model. In this paper, we have 
demonstrated that the data on two others (PG 1115+080, 
Q 2237+030) are consistent with smooth galaxy-like mod- 
els, especially when a realistic treatment of the errors in the 
flux ratios is incorporated. 

We note that the fraction of mass in substructure ex- 
pected in galaxy haloes is very small (< 5 per cent) and it 
occurs overwhelmingly in the outlying portions (e.g., Moore 
ct al. 1999; Ghigna et al. 2000). Substructure evolves as 
it is subjected to tides, impulsive collisions and dynamical 
friction. Tidal disruption becomes important as soon as the 
mean density of the galaxy is equal to the density of the sub- 
structure. Similarly, the dynamical friction timescale scales 
with the square of the galactocentric radius (see e.g., Bin- 
ney & Tremaine 1987). So, both tides and dynamical fric- 
tion are efficient at erasing substructure in the inner parts. 
This effect is referred to as "anti-biasing" by Ghigna et al. 
(2000). So, on dynamical grounds, little substructure is pro- 
jected within the Einstein radius. Taking the example of 
PG 1115+080, the projected Einstein radius is just ~ 3.6 
kpc. The fraction of mass in substructure projected within 
a cylinder of radius ~ 3.6 kpc is clearly much lower than 
the global fraction of 5 per cent, which pertains to the en- 
tire dark halo of total extent ~ 200 kpc. It is crucial that 
lensing calculations do not assume an everywhere uniform 
fraction of substructure in the halo, as this does not take 
into account the "anti-biased" spatial distribution of sub- 
structure and therefore necessarily over-emphasises the im- 
portance of the effects of substructure on flux ratios. Notice 
that an interesting consequence of the spatial distribution is 
that anomalous flux ratios are more likely to occur for lenses 
with larger angular separation, as these have larger Einstein 
radii. 

The real question at issue is the following. Suppose a 
simple lens model (such as a simple isothermal ellipsoid plus 
shear) does not adequately fit the data on positions and 
flux ratios. What can be legitimately deduced? The diffi- 
culty is that there are many modifications of the simple 
model that remove the discrepancy. One of these is sub- 
structure, as pointed out by Dalai & Kochanek and Metcalf 
& Zhao. As shown in this paper, another is higher order mul- 
tipoles in the lensing mass, such as diskiness, boxiness, lop- 
sidedness and barredness (see also Moller, Hewett & Blain 
2003, who make a similar point). It is crucial to develop 
techniques to distinguish between higher order quadrupoles 
on the one hand and substructure on the other. From our 
modelling, we tend to agree that the substructure candi- 
date B 1422+231 originally pointed out by Mao & Schnei- 
der (1998) is strong. Using a novel application of the cusp 
relation, Keeton, Gaudi & Petters (2002) have argued that 
B 2045+265 and RX J 0911+0551 may be two further good 
candidates. 

It is an outstanding problem to predict - for different 
cosmologies - how many quadruplets may have anomalous 
flux ratios. Accurate calculations will become possible only 
when distributions of mass and position of substructure be- 
come available from high resolution simulations. It will be 
interesting to see whether the numbers of lenses for which 
substructure is currently being claimed are compatible with 
cold dark matter theories. 
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APPENDIX A: FULL SOLUTION BY THE 
VARIATION OF PARAMETERS 

Observationally speaking, it is the angular structure of the 
density contours G(9) that is more accessible. It is some- 
times useful to be able to generate the corresponding an- 
gular structure of the lensing potenial F{6). This can be 
done by solving eq (2) using the method of variation of the 
parameters (e.g., Bronshtein & Semendyayev 1998, section 
3.3.1.3.4; Evans & Witt 2001) to give 



F(6) = 



sin(/3fl) 

P 

cos(/36>) 



Ci+ G(#) cos(/?#) di? 



/ 

Jo 



G 2 + / G(tf) sin(/3tf) M 



(Al) 



This establishes the relation between the potential and the 
surface mass density in terms of two constants C\ and G2, 
whose values will be given shortly. 

Therefore, the deflection angle has components 



sin(/3-l)0 



-r^cos^-l)^ 
r^cosCS-l)^ 



Gi + / G(tf)cos(/3tf)dtf 
Jo 

G 2 + [ G(tf) sin(/3i?) dtf 
Jo 



Gi + / G(tf)cos(/3tf)dtf 
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+r p - 1 sm{f3-l)8 



C 2 + 



f 

Jo 



G(tf) sin(/3#) <M 



■ (A2) 



Another expression can be given for the deflection angle, 
namely 



Am 



(A3) 



By substituting the Fourier series into eqs. (A2)-(A3), we 
obtain expressions for the constants 



Ci 

c 2 



kb k 



ao 



(A4) 



This completes the solution written down in eq. (Al). Given 
the angular structure of the density G{9), we can find that 
of the potential F(9), and vice versa. 



APPENDIX B: FOURIER COEFFICIENTS OF 
ELLIPTICAL DISTRIBUTIONS 

Let us consider the Fourier expansion of the elliptical distri- 
bution in the flat rotation curve case (/3 = 1) 



F{&) = A{cos i 6 + q~ sin 



2 • 2 a \-l/2 
cin hi I ' 



The Fourier coefficients are 

i r** 

a k = - F(9)cos(k9)d9 



bk 



1 



2ir 



F(9)sm(k9)d9. 



(Bl) 



(B2) 



(B3) 



However, the only non- vanishing Fourier coefficients are a 2n 
where n is an integer. The lowest-order terms are 



ao = — K 

TV 

4A 

a 2 



7T(l- 9 2 ) 

4A 



(l + q 2 )K-2q 2 E] 



" 4 = Ml _ q2)2 [(^ + l)(q 2 +3)K-8q(l+q 2 )E] 



(Hi 



1 ! [(l + g 2 )(15 + 98g 2 + 15gV (B4) 



f57r(f-g 2 ) 3 

-2g 2 (23+82g 2 + 23g 4 )£] 

° 8 = 1057r(l- g ^4 [ - 32g2(1+g2)(11 + 74g2 + llg4 
+(1 05+f 436g 2 +3062g 4 + 1436g 6 + W5q s )K] 
where K and E denote the following integrals 

/ 2 ^ 



K = / 
Jo 



(1 - (1 -g- 2 )sin 2 tf)V2 



tt/2 



/ = / dti(l - (1 - g~ 2 )sin 2 tf) 1/2 
/o 



(B5 



For elliptical potential models, this gives the Fourier se- 
ries for the angular function F(9). Note that the expres- 
sions become more and more cumbersome with increasing 



order. However, the Fourier series converges (numerically) 
extremely rapidly. 

For elliptical density models, exactly the same coeffi- 
cients hold good but for the Fourier expansion for G(9), the 
shape function in the density. Of course, the Fourier coeffi- 
cients of G(9) are directly related to the Fourier coefficients 
of the potential contour F(9) via eqs. (4) and (5). 



APPENDIX C: CALCULATION OF THE 
ISOPHOTAL DISTORTIONS 

Most giant elliptical galaxies have isophotes that are very 
well approximated by ellipses. The same is of course not true 
for disk or barred or interacting galaxies. Data from numer- 
ical simulations of halo formation are sparse, but suggest 
that dark haloes may also be well-approximated by ellipses 
(e.g., Heyl et al. 1994). So, it is useful to be able to estimate 
the size of the distortions of our fits from pure ellipses. This 
is done by computing Bender et al.'s (1989) coefficients, de- 
fined as 



B 



1 f 27r 
= - / (n(0) - 
n Jo 



r c (9))cos(k8), 



Ok = 



(n(0)-r e (0))sin(fc0), 



(CI) 



(C2) 



where n is the polar equation of the isophote and r e is 
the best fitting ellipse (see also Binney & Merrifield 1998). 
In particular, a positive (or negative) af means that the 
isophotes are disk-like (or boxy). 

From eqn (5), we deduce the equation of the isophote 



n = y + ^[o fc (l~fc 2 ) cos(k9) + b k (l-k 2 ) sin(fc(9)]. (C3) 

fc=2 

Here, we have specialised to the case of a flat rotation curve 
by setting f3 = 1. As shown in Section 3.3, we can transform 
to the coordinate system aligned with the major and minor 
axis and so we can set b 2 = 0. 

Now we must determine the best fitting ellipse, which 
is the one which reproduces the same Fourier coefficients 
ao and a 2 as our isophote. In other words, we solve for the 
amplitude A and axis ratio q in eqn (Bl), 



4A 



ao 



—3a 2 = 



4A 



77(1-02) 



[(l + g 2 )^-2g 2 £]. 



(C4) 



These equations do not have analytic solutions, but the roots 
can be easily obtained numerically. Having solved for A and 
q, our isophote n and the best fitting ellipse r c have Fourier 
series which differ only for coefficients a k and bk with k > 3. 

It is now straightforward to derive formulae for Bender 
ct al.'s coefficients. Tabulated values are available in the 
literature only for the a § and af, so we confine explicit 
formulae to these cases: 



B 
B 

a 4 



-8a 3 , 
-15a 4 - 



4A 



(C5) 

,-[(3g 2 + l)(g 2 +3)^-8g 2 (l+g 2 )iJ] 



37r(l-g 2 )2' 

Numerical values for the lens galaxies in Q 2237+030 
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PG 1115+080 are given in the main text. In both cases, the 
amplitudes of af and af are within the expected ranges, 
confirming that our solutions are realistic. 



APPENDIX D: INCLUSION OF EXTERNAL 
SHEAR 

In this appendix, we give the equations incorporating the 
effects of external shear. The lens equation becomes 

£ = x+^x+yzy- — , ri = y + yzx-^y-—. (Dl) 
The vector of observables becomes 
d = ^(l+7i)rf cos0* +72r-£sin0f, . . . , 
(I-71 )r* sin 6^ + 72 r* cos 0*, ... , 

(/«-l)r m (l- 7 i-72 2 ),...) T , (D2) 

where I runs from 1 to n for a lens system with n images. 
The matrix C retains the same form as given in eq (22), but 
the coefficients ~fk£ and Ski become 

lk{0t) = (l-fe 2 )[/ w r,cos(fce 1 )H/(^ 1 )-r 1 cos(fc^)W(^)] 
5 k {0e) = (l-k 2 )[f le r e sm(ke 1 )W(6 1 )-r 1 sm(ke e )W(9 e )], 

with 

W{6) = I+71 cos 20* +72 sin 20*. (D3) 

The vector of unknowns is unchanged from eq (21). For a 
given external shear, the problem remains linear and can be 
solved by SVD. 
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